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Abstract 

Taylor impacts tests were originally devised to determine the dynamic yield strength of materials at moderate 
strain rates. More recently, such tests have been used extensively to validate numerical codes for the simulation 
of plastic deformation. In this work, we use the material point method to simulate a number of Taylor impact tests 
to compare different Johnson-Cook, Mechanical Threshold Stress, and Steinberg-Guinan-Cochran plasticity 
models and the vob Mises and Gurson-Tvergaard-Needleman yield conditions. In addition to room temperature 
Taylor tests, high temperature tests have been performed and compared with experimental data. 

1 INTRODUCTION 

The Taylor impact test (Taylor [ Q) was originally devised as a means of determining the dynamic yield strength of 
solids. The test involves the impact of a fiat-nosed cylindrical projectile on a hard target at normal incidence. Taylor 
provided an analytical solutions for the dynamic yield strength of the material of the projectile based on the length 
of the elastic region and the radius of the region of permanent set. As described by Whifhn O, that use of the test 
was limited to relatively small deformations obtained from low velocity impacts. Though the Taylor impact test 
continues to be used to determine yield strengths of materials at high strain rates, the test is limited to peak strains 
of around 0.6 at the center of the specimen (Johnson and Holmquist J3J). For higher strains and strain rates, the 
Taylor test is currently used more as a means of validating plasticity models in numerical codes for the simulation 
of high rate phenomena such as impact and explosive deformation as suggested by Zerilli and Armstrong 0. 

In this paper, we describe our experience in validating the plasticity models in a parallel, multiphysics 
code that uses the material point method (Sulsky et al. 016]]) using Taylor impact tests for various strain rates and 
temperatures. A number of metrics are used to compare simulations and experiments and suggesstions are made 
regarding the use of Taylor impacts tests for the validation of the plasticity portion of such codes. 

The organization of this paper is as follows. Section [2] provides the background for the current study and 
describes the mutiphysics code Uintah, the material point method, and the stress update algorithm, and various 
plasticity models and yield conditions. A few validation metrics are identified and their significance is discussed 
in Section [3] Comparisons between experimental data and simulations of Taylor impact tests using the validation 
metrics are described in Section [4] Finally, conclusions and suggestions are presented in Section [5] 

2 BACKGROUND 

The goal of this work is to present some results and insights we have obatined during the process of validation 
of plasticity models used in the simulation of the deformation and failure of a steel container that expands under 
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Figure 1: Simulation of exploding cylinder. 



the effect of gases produced by an explosively reacting high energy material (PBX 9501) contained inside. The 
entire process is simulated using the massively parallel, Common Component Architecture Q based, Uintah 
Computational Framework (UCF) (HJ . 

The high energy material reacts at temperatures of 450 K and above, This elevated temperature is achieved 
through external heating of the steel container. Experiments conducted at the University of Utah have shown that 
failure of the container can be due to ductile fracture associated with void coalescence and adiabatic shear bands. 
If shear bands dominate the steel container fragments, otherwise a few large cracks propagate along the cylinder 
and pop it open. Figure [T] shows the result of a simulation of a coupled fire-container-explosion using the Uintah. 

The dynamics of the solid materials - steel and PBX 9501 - is modeled using the Lagrangian Material 
Point Method (MPM) [5]. Gases are generated from solid PBX 9501 using a burn model [9]. Gas-solid interaction 
is accomplished using an Implicit Continuous Eulerian (ICE) multi-material hydrodynamic code |[T0l . A single 
computational grid is used for all the materials. 

The constitutive response of PBX 9501 is modeled using ViscoSCRAM [11], which is a five element 
generalized Maxwell model for the viscoelastic response coupled with statistical crack mechanics. Solid PBX 
9501 is progressively converted into a gas with an appropriate equation of state. The temperature and pressure in 
the gas increase rapidly as the reaction continues. As a result, the steel container is pressurized, undergoes plastic 
deformation, and finally fragments. 

The main issues regarding the constitutive modeling of the steel container are the selection of appropriate 
models for nonlinear elasticity, plasticity, damage, loss of material stability, and failure. The numerical simulation 
of the steel container involves the choice of appropriate algorithms for the integration of balance laws and consti- 
tutive equations, as well as the methodology for fracture simulation. Models and simulation methods for the steel 
container are required to be temperature sensitive and valid for large distortions, large rotations, and a range of 
strain rates (quasistatic at the beginning of the simulation to approximately 10 6 s _1 at fracture). 

The approach chosen for the present work is to use hypoelastic-plastic constitutive models that assume an 
additive decomposition of the rate of deformation tensor into elastic and plastic parts. Hypoelastic materials are 
known not to conserve energy in a loading-unloading cycle unless a very small time step is used. However, the 
choice of this model is justified under the assumption that elastic strains are expected to be small for the problem 
under consideration and unlikely to affect the computation significantly. 

Two plasticity models for flow stress are considered along with a two different yield conditions. Explicit 
fracture simulation is computationally expensive and prohibitive in the large simulations under consideration. The 
choice, therefore, has been to use damage models and stability criteria for the prediction of failure (at material 
points) and particle erosion for the simulation of fracture propagation. 
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2.1 The Material Point Method 

The Material Point Method (MPM) [51 is a particle method for structural mechanics simulations. In this method, 
the state variables of the material are described on Lagrangian particles or "material points". In addition, a regular, 
structured Eulerian grid is used as a computational scratch pad to compute spatial gradients and to solve the gov- 
erning conservation equations. An explicit time-stepping version of the Material Point Method has been used in 
the simulations presented in this paper. The MPM algorithm is summarized below 

It is assumed that an particle state at the beginning of a time step is known. The mass (m), external force 
(f ext ), and velocity (v) of the particles are interpolated to the grid using the relations 

m 9 = Yl S 9P m P ' V f = (V m s) S 9P m P V P ' ^ = Y S 9P ^ W 
P P P 

where the subscript (g) indicates a quantity at a grid node and a subscript (p) indicates a quantity on a particle. The 
symbol ^ indicates a summation over all particles. The quantity (S gp ) is the interpolation function of node (g) 
evaluated at the position of particle (p). Details of the interpolants used can be found elsewhere lfl2l . 

Next, the velocity gradient at each particle is computed using the grid velocities using the relation 

Vv p = £ G gp v g (2) 

9 

where G gp is the gradient of the shape function of node (g) evaluated at the position of particle (p). The velocity 
gradient at each particle is used to determine the Cauchy stress (er p ) at the particle using a stress update algorithm. 
The internal force at the grid nodes (f™ 1 ) is calculated from the divergence of the stress using 

f <? nt = ^9P a pVp (3) 

p 

where V p is the particle volume. 

The equation for the conservation of linear momentum is next solved on the grid. This equation can be cast 
in the form 

m, a g = ff - fj* (4) 

where & g is the acceleration vector at grid node (g). 

The velocity vector at node (g) is updated using an explicit (forward Euler) time integration, and the particle 
velocity and position are then updated using grid quantities. The relevant equations are 

v 9 (t + At) = v 9 (t) + a 9 At (5) 
v p (t + At)=v p (t) + ^2S gp a g At; x p (t + At) = x p (t) + ^ S gp v g At (6) 

9 9 

The above sequence of steps is repeated for each time step. The above algorithm leads to particularly 
simple mechanisms for handling contact. Details of these contact algorithms can be found elsewhere [13 ]. 



2.2 Plasticity and Failure Simulation 

A hypoelastic-plastic, semi-implicit approach lfT4l has been used for the stress update in the simulations presented 
in this paper. An additive decomposition of the rate of deformation tensor into elastic and plastic parts has been 
assumed. One advantage of this approach is that it can be used for both low and high strain rates. Another advantage 
is that many strain-rate and temperature-dependent plasticity and damage models are based on the assumption of 
additive decomposition of strain rates, making their implementation straightforward. 

The stress update is performed in a co-rotational frame which is equivalent to using the Green-Naghdi 
objective stress rate. An incremental update of the rotation tensor is used instead of a direct polar decomposition 
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of the deformation gradient. The accuracy of model is good if elastic strains are small compared to plastic strains 
and the material is not unloaded. It is also assumed that the stress tensor can be divided into a volumetric and a 
deviatoric component. The plasticity model is used to update only the deviatoric component of stress assuming 
isochoric behavior. The hydrostatic component of stress is updated using a solid equation of state. 

Since the material in the container may unload locally after fracture, the hypoelastic -plastic stress update 
may not work accurately under certain circumstances. An improvement would be to use a hyperelastic-plastic 
stress update algorithm. Also, the plasticity models are temperature dependent. Hence there is the issue of severe 
mesh dependence due to change of the governing equations from hyperbolic to elliptic in the softening regime 
IPT51 IT6l IT71 . Viscoplastic stress update models or nonlocal/gradient plasticity models IT8l |T9ll can be used to 
eliminate some of these effects and are currently under investigation. 

A particle is tagged as "failed" when its temperature is greater than the melting point of the material at the 
applied pressure. An additional condition for failure is when the porosity of a particle increases beyond a critical 
limit. A final condition for failure is when a bifurcation condition such as the Drucker stability postulate is satisfied. 
Upon failure, a particle is either removed from the computation by setting the stress to zero or is converted into 
a material with a different velocity field which interacts with the remaining particles via contact. Either approach 
leads to the simulation of a newly created surface. 

In the parallel implementation of the stress update algorithm, sockets have been added to allow for the 
incorporation of a variety of plasticity, damage, yield, and bifurcation models without requiring any change in the 
stress update code. The algorithm is shown in Algorithm[T] The equation of state, plasticity model, yield condition, 
damage model, and the stability criterion are all polymorphic objects created using a factory idiom in C++ [20]. 



2.3 Models 



The stress in the solid is partitioned into a volumetric part and a deviatoric part. Only the deviatoric part of stress 
is used in the plasticity calculations assuming isoschoric plastic behavior. 

The hydrostatic pressure (p) is calculated either using the bulk modulus (K) and shear modulus (/x) or from 
a temperature-corrected Mie-Gruneisen equation of state of the form |[T4l 



P 



PoClC [1 + (1 



[i - (s a - i)cr + r c p r 



C = (p/po - 1) 



(7) 



where Co is the bulk speed of sound, po is the initial density, p is the current density, C p is the specific heat 
at constant volume, T is the temperature, To is the Gruneisen's gamma at reference state, and S a is the linear 
Hugoniot slope coefficient. 

Depending on the plasticity model being used, the pressure and temperature-dependent shear modulus (p) 
and the pressure-dependent melt temperature (T m ) are calculated using the relations 11211 



P = Po 



l + A 



P 

1/3 



T m = T m0 exp 



2a [ 1 



B(T - 300) 



2(r -a-l/3) 



(8) 
(9) 



where, po is the shear modulus at the reference state(T = 300 K, p = 0, e p = 0), e p is the plastic strain, rj = p/po is 
the compression, A = (1/ po)(dp/dp), B = (1/ po)(dp/dT), T m o is the melt temperature at p = po, and a is the 
coefficient of the first order volume correction to Gruneisen's gamma. 

We have explored two temperature and strain rate dependent plasticity models - the Johnson-Cook plasticity 
model (221 and the Mechanical Threshold Stress (MTS) plasticity model lEUHU. The flow stress (07) from the 
Johnson-Cook model is given by 



a f = [A + B(e p ) n ][l + CHe;)][l-(T* 



€p0 



{T-T r 



(10) 
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Data: Persistent:Initial moduli, temperature, porosity, scalar damage, equation of state, plasticity model, 

yield condition, stability criterion, damage model 
Temporary:Particle state at time t 
Result: Particle state at time t + At 
for all the patches in the domain do 

Read the particle data and initialize updated data storage; 
for all the particles in the patch do 

Compute the velocity gradient, the rate of deformation tensor and the spin tensor; 
Compute the updated left stretch tensor, rotation tensor, and deformation gradient; 
Rotate the input Cauchy stress and the rate of deformation tensor to the material configuration; 
Compute the current shear modulus and melting temperature; 

Compute the pressure using the equation of state, update the hydrostatic stress, and compute the trial 
deviatoric stress; 

Compute the flow stress using the plasticity model; 
Evaluate the yield function; 
if particle is elastic then 

Rotate the stress back to laboratory coordinates; 
Update the particle state; 

se 

Find derivatives of the yield function; 
Do radial return adjustment of deviatoric stress; 

Compute updated porosity, scalar damage, and temperature increase due to plastic work; 
Compute elastic-plastic tangent modulus and evaluate stability condition; 
Rotate the stress back to laboratory coordinates; 
Update the particle state; 

if Temperature > Melt Temperature or Porosity > Critical Porosity or Unstable then 

j Tag particle as failed; 
end 

end 

end 
end 

Convert failed particles into a material with a different velocity field; 

Algorithm 1: Stress Update Algorithm 



where e p o is a user defined plastic strain rate, A, B, C, n, m are material constants, T r is the room temperature, and 
T m is the melt temperature. 

The flow stress for the MTS model is given by 



where 



■ M o ~ . M Q * 

Pf = 0~a-\ <->iCi "I ^eCe 

Mo Mo 



(11) 



D 



M = Mo 



Si 



exp(f; 



kT epj 
goilib 3 e 



l/qi 



1/pi 



kT 



l/qe 
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e = e [l- F(X)} + 9 IV F(X) ; e = a + ai In e + a 2 Ve - a 3 T 
X = ^ ; = tanh(aX) ; ln(«W<7 es0 ) = ( -^—) In [ — 

Oes \lur goes J ytesO^ 



and a a is the athermal component of mechanical threshold stress, fiQ is the shear modulus at K, D, Tq are em- 
pirical constants, &i represents the stress due to intrinsic barriers to thermally activated dislocation motion and 
dislocation-dislocation interactions, a e represents the stress due to microstructural evolution with increasing de- 
formation, k is the Boltzmann constant, b is the length of the Burger's vector, <7o[i,e] are the normalized activation 
energies, £o[j,e] are constant strain rates, q^ e j , p[ i>e ] are constants, 9q is the hardening due to dislocation accumu- 
lation, ao, a\, a,2, as, OiVi ol are constants, a es is the stress at zero strain hardening rate, a es o is the saturation 
threshold stress for deformation at K, go es is a constant, and e es o is the maximum strain rate. 

We have decided to focus on ductile failure of the steel container. Accordingly, two yield criteria have been 
explored - the von Mises condition and the Gurson-Tvergaard-Needleman (GTN) yield condition ||25l |26ll which 
depends on porosity. An associated flow rule is used to determine the plastic rate parameter in either case. The von 
Mises yield condition is given by 



^) 2 -l = °; 



$ = ( — ) -1 = 0; a eq = \l-a d :a d (12) 



where a eq is the von Mises equivalent stress, a d is the deviatoric part of the Cauchy stress, and a? is the flow stress. 
The GTN yield condition can be written as 

d> = (^j 2 + 2q % U cosh (<&^) " (1 + = (13) 
where q\ , q%, qs are material constants and /* is the porosity (damage) function given by 

/* = r for ; - /c ' d4) 

' ' ■ - f c ) for f>f c 

where k is a constant and / is the porosity (void volume fraction). The flow stress in the matrix material is 
computed using either of the two plasticity models discussed earlier. Note that the flow stress in the matrix material 
also remains on the undamaged matrix yield surface and uses an associated flow rule. 

The evolution of porosity is calculated as the sum of the rate of growth and the rate of nucleation [27]. The 
rate of growth of porosity and the void nucleation rate are given by the following equations ESI 




/ — /nucl + /grow (15) 

/grow = (1 - /)Tr(D p ) (16) 



; fn 
/nucl = 7==7 ex P 



1 (gp ~ £n) 2 

2^T 



(17) 



where D p is the rate of plastic deformation tensor, f n is the volume fraction of void nucleating particles , e n is the 
mean of the distribution of nucleation strains, and s n is the standard deviation of the distribution. 

Part of the plastic work done is converted into heat and used to update the temperature of a particle. The 
increase in temperature (AT) due to an increment in plastic strain (Ae p ) is given by the equation |29l 



(18) 
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where x is the Taylor-Quinney coefficient, and C p is the specific heat. A special equation for the dependence of C p 
upon temperature is also used for steel QUI . 



C p = 10 3 (0.09278 + 7.454 x 1(T 4 T + 12404. 0/T 2 ) (19) 

Under normal conditions, the heat generated at a material point is conducted away at the end of a time step 
using the heat equation. If special adiabatic conditions apply (such as in impact problems), the heat is accumulated 
at a material point and is not conducted to the surrounding particles. This localized heating can be used to simulate 
adiabatic shear band formation. 

After the stress state has been determined on the basis of the yield condition and the associated flow rule, a 
scalar damage state in each material point can be calculated using either of two damage models - the Johnson-Cook 
model Oil or the Hancock-MacKenzie model [32]. While the Johnson-Cook model has an explicit dependence 
on temperature, the Hancock-McKenzie model depends on the temperature implicitly, via the stress state. Both 
models depend on the strain rate to determine the value of the scalar damage parameter. 

The damage evolution rule for the Johnson-Cook damage model can be written as 



D 3 

D 1 + D 2 exp —a* 



[1 + L> 4 ln(e p *)] [1 + D 5 T*] ; a* = ; (20) 



eq 



where D is the damage variable which has a value of for virgin material and a value of 1 at fracture, e{, is the 
fracture strain, Z?i, D2, D%, D±, D§ are constants, cr is the Cauchy stress, and T* is the scaled temperature as in the 
Johnson-Cook plasticity model. 

The Hancock-MacKenzie damage evolution rule can be written as 

b / = 1-65 

e f ' ^ exp(1.5a*) 

The determination of whether a particle has failed can be made on the basis of either or all of the following 
conditions: 

• The particle temperature exceeds the melting temperature. 

• The TEPLA-F fracture condition [33] is satisfied. This condition can be written as 

(f/fcf + (e P /el) 2 = 1 (22) 

where / is the current porosity, f c is the maximum allowable porosity, e p is the current plastic strain, and Cp 
is the plastic strain at fracture. 

• An alternative to ad-hoc damage criteria is to use the concept of bifurcation to determine whether a particle 
has failed or not. Two stability criteria have been explored in this paper - the Drucker stability postulate Q4l 
and the loss of hyperbolicity criterion (using the determinant of the acoustic tensor) 05ll36l . 

The simplest criterion that can be used is the Drucker stability postulate 04l which states that time rate of 
change of the rate of work done by a material cannot be negative. Therefore, the material is assumed to become 
unstable (and a particle fails) when 

&:V P <0 (23) 

Another stability criterion that is less restrictive is the acoustic tensor criterion which states that the material 
loses stability if the determinant of the acoustic tensor changes sign ||35l[36). Determination of the acoustic tensor 
requires a search for a normal vector around the material point and is therefore computationally expensive. A 
simplification of this criterion is a check which assumes that the direction of instability lies in the plane of the 
maximum and minimum principal stress 071 . In this approach, we assume that the strain is localized in a band 
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with normal n, and the magnitude of the velocity difference across the band is g. Then the bifurcation condition 
leads to the relation 

Rijgj = ; Rij = M ikj in k ni + M ilkj n k ni - a ik rijn k (24) 

where are the components of the co-rotational tangent modulus tensor and Oij are the components of the 

co-rotational stress tensor. If det(-Ry) < 0, then gj can be arbitrary and there is a possibility of strain localization. 
If this condition for loss of hyperbolicity is met, then a particle deforms in an unstable manner and failure can be 
assumed to have occurred at that particle. 

3 VALIDATION METRICS 

The attractiveness of the Taylor impact test arises because of the simplicity and inexpensiveness of the test. A 
flat-ended cylinder is fired on a target at a large enough velocity and the final deformed shape is measured. The 
drawback of this test is that intermediate states of the cylinder are difficult to measure and hence are generally 
not. The validation metrics that we consider in this paper are based on the final shape of the cylinder though other 
metrics may be considered if measurements of these are made during the course of an impact test. We note that the 
Taylor test could also be used to validate simulations of dynamic fracture though we do not address that issue in 
this paper. 

There is a large literature on the systematic verification and validation of computational codes (see Oberkampf 
et al. ll38ll . Babuska and Oden [39] and references therein). It has been suggested that validation metrics be devel- 
oped that can be used to compare experimental data and simulation results. The metrics discussed in this paper are 
intended to be a step in that direction but they are not intended to be complete or comprehensive. 

The most common metric used in the literature is the "calibrated eyeball" approach or "view-graph norm" 
(Oberkampf et al. ||38ll ) where a plot of the simulated deformed configuration is superimposed on the experimental 
data and a subjective judgement of accuracy is made. We believe that there is value to this approach and present 
all our data in this form. However, we also believe that more quantitative descriptions of the difference between 
experiment and simulations can be obtained and present comparisons using other metrics. 

Metrics, sensitivity studies, and determination of experimental variability are essential. Some quantities of 
interest are: 



Metrics 


(a) 


Regression between profiles 


(b) 


Length change 


(c) 


Diameter change 


(d) 


Volume change 


(e) 


Middle bulge difference 


(f) 


Length of elastic zone 


(g) 


plastic strain 


(h) 


temperature 


(i) 


time of impact 


0) 


energy conversion at impact 



2. Sensitivity studies 

(a) mesh size (quantify discretization errors) 

(b) plasticity model parameters 

(c) plasticity model 
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(d) impact velocity 

(e) temperature 

(f) length and diameter 

3. Variability in experimental data 

(a) material 

(b) geometry 

(c) velocity 

(d) temperature 

(e) measurement error 



4 Taylor impact simulations 



In this section, we compare the final deformed shapes from simulations of Taylor impact tests with experimentally 
obtained data. In cases where images or profiles of the deformed shapes of the cylinders were available, these 
were digitized using a scanner and then imported into XFig (Sutanthavibul et al. Il40l ). The scanned images were 
overlaid with manually digitized lines that were drawn as accurately as possible after expanding the images to a 
resultion of 1024x 1260. The digitized curves were then rotated so that the axes were aligned with the grid. The 
XFig coordinates were then scaled to length units using cues from the digitized images and their axes (if any were 
provided). Some small errors (l%-2%) are expected in this procedure. However, the overall profiles of the cylinders 
are captured accurately in most cases. 

The simulations were run for 150 /us - 200 fis depending on the problem. The simulation times were chosen 
such that the cylinders bounced off the anvil and moved away for at least 20 fj£. It was observed at beyond this 
time, the deformed shape of the cylinder reamined constant and all elastic strains and rotations had been recovered. 

In the paper by Carrington and Gayler [4 1 ] a highly deformed mild steel specimen has been shown (plate 
1, figure 3). To determine if MPM could be used to simulate such large deformations, we ran a Taylor impact test 
on the problem geometry using the Johnson-Cook plasticity model for 4340 steel. The final deformed shape from 
Carrington's paper is compared with our predicted shape in Figure [2j The initial velocity is Vq = 2140ft/ s = 
652.3m/s, the initial diameter of the cylinder is Dq = 0.5m = 12.7mm, and the initial length of the cylinder is 
L = O.mOin = 25.37mm. 




(a) Actual profile (Carrington and Gayler fiTTl ). 



Expt. 
MPM 




-0.5 0.5 
cm 



(b) Computed vs. actual profile. 



Figure 2: Comparison of experimental vs. computed shapes. L = 25.37 mm , D = 12.7 mm, Vq = 652.3 m/s. 



In our simulation, the 4340 steel cylinder was impacted against a stiff anvil using frictional contact. The 
4340 steel flows much more readily than the mild steel used by Carrington and Gayler [41]. The experiment also 
shows that the tips of the "mushroom" have broken off. We did not simulate any fracture and hence we do not see 
that effect. However, the overall shape of the deformed specimen suggests that our simulations can provide good 
qualitative descriptions of large deformations. To quantify how well our simulations fit experimental data, we ran 
a series of Taylor impact tests on various materials and compared them against experimental data. Some of those 
results are presented in this report. 
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Table 1: Initial data for copper simulations. OFHC = oxygen free high conductivity. ETP = electrolytic tough pitch. 



Case 


Material 


Initial 


Initial 


Initial 


Initial 


Source 






Length 


Diameter 


Velocity 


Temperature 








(L mm) 


(D mm) 


(Vo m/s) 


(To K) 




Cu-A 


OFHC Cu 


23.47 


7.62 


210 


298 


Wilkins and Guinan [42 ] 


Cu-B 


OFHC Cu 


25.4 


7.62 


130 


298 


Johnson and Cook E21 


Cu-C 


OFHC Cu 


25.4 


7.62 


146 


298 


Johnson and Cook ll22l 


Cu-D 


OFHC Cu 


25.4 


7.62 


190 


298 


Johnson and Cook [22] 


Cu-E 


E1P Cu 


31) 


6.00 


111 


one 
295 


uust 1143 II 


Cu-F 


ETPCu 


30 


6.00 


188 


718 


Gust El 


Cu-G 


ETPCu 


30 


6.00 


211 


727 


Gust El 


Cu-H 


ETP Cu 


30 


6.00 


178 


1235 


Gust El 


Cu-I 


Annealed Cu. 










Zocher et al. HH 


Cu-J 


With porosity 










Addessio et al. EH 



4.1 Taylor impact tests on copper 

In this section we present the results from Taylor tests on copper specimens for different initial temperatures and 
impact velocities. Table[T]shows the initial dimensions, velocity, and temperature of the specimens (along with the 
type of copper used and the source of the data) that we have simulated and compared with experimental data. 



4.1.1 Room temperature impact of copper 

Comparisons between the computed and experimental profiles of annealed copper specimen Cu-I are shown in 
Figure [3] The MTS model predicts the final length quite accurately (at this is true for other room temperature 
simulations of copper). The profile shape is also computed accurately. The Johnson-Cook model overestimate the 
final length. However, the difference is small and may be attributed to material variability. 



4 
3.5 

3 
2.5 

2 
1.5 

1 

0.5 



-1 1 

(a) Johnson-Cook. 



4 

3.5 
3 

2.5 
2 

1.5 
1 

0.5 




1 

(b) Mechanical Threshold Stress. 



Figure 3: Comparison of experimental and computed shapes of annealed copper cylinder Cu-I using the Johnson- 
Cook and Mechanical Threshold Stress plasticity models. The axes are shown in cm units. 

Simulations of impact case Cu-I with increasing mesh refinement are shown in Figure |4| The number of MPM 
particles is doubled with each refinement. We observe that the solution does not change much as we refine the 
mesh. However, this is true only at low temperatures and moderate impact velocities. Significant mesh dependence 
is observed at high temperatures where softening becomes dominant as wee will see in our calculations with 6061- 
T6 aluminum. 
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Figure 4: Comparison of experimental and computed shapes of 606 1T6 aluminum cylinders using the Johnson- 
Cook (JC) with increasing mesh refinement. The axes are in cm. 



4.1.2 High temperature impact of copper 

At higher temperatures, the response of the three plasticity models is quite different. Comparisons between the 
computed and experimental profiles of ETP copper specimen Cu-F are shown in Figure [5j a), (b), and (c). Those 
for specimen Cu-G are shown in Figure [5jd), (e), and (f). If frictional contact at the impact surface is simulated, 
the final shapes of the specimens Cu-F and Cu-G are as shown in Figure|5Jg), (h), (i), (j), (k), and (1). 

Notice that though both specimens are nominally at the same temperature and has almost identical impact 
velocities, the final profile is quite different even though the final lengths are nearly identical. It is likely that most 
of the difference is due the initial conditions with a small contribution from material variability. This conjecture is 
partially supported by the fact that the profiles predicted by the Johnson-Cook model match the experiments quite 
well. 

We observe that the Johnson-Cook and Steinberg-Guinan models perform well for specimen Cu-F when friction 
is not included in the calculation. In the presence of frictional contact, the predicted profiles deviate significantly 
from the experimental profiles in the mushroom region. This indicates that there is a possibility of inaccurate 
contact force calculation when friction is included. 

From specimen Cu-G, the slightly higher impact velocity leads to an underestimation of the final length by the 
Johnson-Cook model, even though the mushroom region is predicted accurately. The MTS model overestimates 
the length and underestimates the mushroom diameter while the Steinberg-Guinan model predicts the final length 
best but fails to predict the mushroom shape. Once again, frictional contact appears to reduce the accuracy of the 
prediction. 



4.1.3 Comparisons with FEM 

To determine how our MPM simulations compare with FEM simulations we have run two high temperature ETP 
copper impact tests using LS-DYNA (with the coupled structural-thermal option). Figure [6] shows the final de- 
formed shapes for the two cases from the MPM and FEM simulations using Johnson-Cook plasticity. In this case 
frictional contact has been used. 

The FEM simulations consistently overestimate the final length of the specimen though the mushroom diameter 
is more accurately predicted by FEM. For the case where no contact friction is applied, MPM predictions are 
consistently better than FEM predictions. 
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(a)JC(Cu-F). (b) MTS (Cu-F). (c) SCG (Cu-F). (d) JC (Cu-G). (e) MTS (Cu-G). (f) SCG (Cu-G). 








(g) JC (Cu-F) 
with friction. 



(h) MTS (Cu-F) 
with friction. 



(i) SCG (Cu-F) 
with friction. 



(j) JC (Cu-G) 
with friction. 



(k) MTS (Cu-G) (1) SCG (Cu-G) 
with friction. with friction. 



Figure 5: Comparison of experimental and computed shapes of ETP copper cylinders using the Johnson-Cook 
(JC), Mechanical Threshold Stress (MTS), and Steinberg-Cochran-Guinan (SCG) plasticity models. Specimen C-F 
has an initial temperature of 718 K and Cu-G is initially at 727 K. The initial velocities are 188 m/s and 211 m/s, 
respectively. The axes are shown in cm units. 




Figure 6: Comparison of experimental and computed shapes of ETP copper cylinders using MPM and FEM. The 
axes are in cm. 



C-SAFE-CD-IR-04-004 



13 



4.2 Taylor impact tests on 6061-T6 aluminum alloy 

In this section we present the results from Taylor tests on 6061-T6 aluminum specimens for different initial temper- 
atures and impact velocities. We have chosen to study this material as it is a well characterized face centered cubic 
material that has been utilized by Chhabildas et al. ||45l for the validation of high velocity impacts that formed the 
basis of the second stage of our validation simulations. Table [2] shows the initial dimensions, velocity, and temper- 
ature of the specimens (along with the type of copper used and the source of the data) that we have simulated and 
compared with experimental data. 

Table 2: Initial data for 6061-T6 aluminum simulations. 



Case 


Material 


Initial 


Initial 


Initial 


Initial 


Source 








Length 


Diameter 


Velocity 


Temperature 










(L Q mm) 


(D mm) 


(Vb m/s) 


(7b K) 






Al-A 


6061-T6 Al 


23.47 


7.62 


373 


298 


Wilkins and Guinan 


El 


Al-B 


6061-T6 Al 


23.47 


7.62 


603 


298 


Wilkins and Guinan 


an 


Al-C 


6061-T6 Al 


46.94 


7.62 


275 


298 


Wilkins and Guinan 


ma 


Al-D 


6061-T6 Al 


46.94 


7.62 


484 


298 


Wilkins and Guinan 




Al-E 


6061-T6 Al 


30 


6.00 


200 


295 


Gust 1431 




Al-F 


6061-T6 Al 


30 


6.00 


358 


295 


GustHa 




Al-G 


6061-T6 Al 


30 


6.00 


194 


635 


Gust ma 




Al-H 


6061-T6 Al 


30 


6.00 


354 


655 


Gustga 




Al-I 


6061-T6 Al 










Addessio et al. O 





4.2.1 Room temperature impact: 6061-T6 Al 

Comparisons between the computed and experimental profiles of 606 1T6 aluminum alloy specimen Al-A are 
shown in Figure |7ja), (b), and (c). Those for specimen Al-C are shown in Figure |7Jd), (e), and (f). If frictional 
contact at the impact surface is simulated, the final shapes of the specimens Al-A and Al-C are as shown in Figure[7] 

(g), (h), (i), (j), (k), and (1). 

We note that all three models predict essentially identical profiles. The higher velocity impact of the shorter 
specimen Al-A is best predicted by the MTS model as far as final length is concerned. The mushroom width is 
predicted better when some friction is included at the anvil-specimen interface. There is a noticeable amount of 
curvature under frictional contact. We believe that this partly due to the contact algorithm that has been used. 

The longer specimens have lower impact velocities. However, all three models predict a final length that is 
shorter than that observed in experiment. We believe that is discrepancy is due to material variability. Note the 
accuracy with which the profiles are predicted and the noticeably lower curvature of the mushroom under frictional 
contact compared to specimen Al-A. 

4.2.2 High temperature impact: 6061-T6 Al 

At higher temperatures, the response of the three plasticity models is quite different. Comparisons between the 
computed and experimental profiles of 606 1T6 aluminum alloy specimens have been performed under conditions 
of frictional contact. The final shapes of the specimens Al-G and Al-H are as shown in Figure[8] If failure simulation 
is included, the profiles are as shown in Figures [8jg), (h), (i), (j), (k), and (1). 

For the lower impact velocity of specimen Al-G, the Johnson-Cook model performs the best at predicting both 
the final length and the mushroom diameter. Both the MTS and SGC models overestimate the final length and 
underestimate the mushroom diameter. The MTS model fares slightly worse than the SCG model. However, the 
differences are small enough that they can be attributed to material variability. Including erosion effects in the 
simulation does not affect the result significantly. 
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(a)JC(Al-A). (b) MTS (Al-A). (c) SCG (Al-A). (d) JC (Al-C). (e) MTS (Al-C). (f) SCG (Al-C). 



(g) JC (Al-A) 
Friction. 



(h) MTS (Al-A) (i) SCG (Al-A) 
Friction. Friction. 



<j)JC (Al-C) 
Friction. 



(k) MTS (Al-C) (1) SCG (Al-C) 
Friction. Friction. 



Figure 7: Comparison of experimental and computed shapes of 606 1T6 aluminum cylinders using the Johnson- 
Cook (JC), Mechanical Threshold Stress (MTS), and Steinberg-Cochran-Guinan (SCG) plasticity models. The 
figure in the top row are from simulations without friction while those in the bottom row are with friction. The axes 
are shown in cm units. 



At the higher impact velocity represented by specimen Al-H, all models fail to predict the final length ac- 
curately. The Johnson-Cook model comes closest but overestimates the length and has an excessively deformed 
mushroom region. The MTS and SCG models have more reasonably shaped mushroom regions but fail to predict 
the final length by almost 100%. The SCG model is slightly better than the MTS model. 

It is possible that the discrepancy that we observe for specimen Al-H is due to inadequate discretization. Simu- 
lations of impact specimen Al-H with increasing mesh refinement are shown in Figure |9j The number of grid cells 
in the plane of the specimen profile has been doubled with each refinement. 

If we examine the profiles shown in Figure |9ja), we observe that the cylinder does appear to shorten with 
increasing refinement. However, there is unphysical curling of the end of the specimen. On the other hand, if 
we eliminate friction from the calculation, the mushroom appears to increase with increased refinement while the 
length decreases. This indicates that there is some amount of mesh dependence of the solution that is probably due 
to the softening behavior of the material. 



4.2.3 Comparisons with FEM 

To determine how our MPM simulations compare with FEM simulations we have run two high temperature alu- 
minum impact tests using LS-DYNA (with the coupled structural-thermal option). Figure 10 shows the final de- 
formed shapes for the two cases from the MPM and FEM simulations using Johnson-Cook plasticity. The FEM 
simulations consistently overestimate the final length and underestimate the mushroom diameter at high tempera- 
tures. 
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(a)JC(Al-G). (b) MTS (Al-G). (c) SCG (Al-G). (d) JC (Al-H). (e) MTS (Al-H). (f) SCG (Al-H). 






\ 



5 u 

q| 



(g) JC (Al-G) 
with erosion. 



(h) MTS (Al-G) 
with erosion. 



(i) SCG (Al-G) (j) JC (Al-H) with (k) MTS (Al-H) 
with erosion. erosion. with erosion. 



(1) SCG (Al-H) 
with erosion. 



Figure 8: Comparison of experimental and computed shapes of 606 1T6 aluminum cylinders using the Johnson- 
Cook (JC), Mechanical Threshold Stress (MTS), and Steinberg-Cochran-Guinan (SCG) plasticity models. Speci- 
mens Al-G and Al-H are both initially at 635 K. Al-G has an impact velocity of 194 m/s while Al-H impacts at 354 
m/s. The axes are shown in cm units. 
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Figure 9: Comparison of experimental and computed shapes of 606 1T6 aluminum cylinders (Al-H) using the 
Johnson-Cook (JC) with increasing mesh refinement. The axes are in cm. 
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Figure 10: Comparison of experimental and computed shapes of 606 1T6 aluminum cylinders using MPM and 
FEM. The axes are in cm. 



C-SAFE-CD-IR-04-004 



17 



4.3 Taylor impact tests on 4340 steel 

In this section we present the results from Taylor tests on 4340 steel specimens for different initial temperatures 
and impact velocities. Table[3]shows the initial dimensions, velocity, and temperature of the specimens (along with 
the type of copper used and the source of the data) that we have simulated and compared with experimental data. 
Note that only a few representative results are shown in this report. 

Table 3: Initial data for 4340 steel simulations. 



Case 


Hardness 


Initial 
Length 
(L mm) 


Initial 
Diameter 
(D mm) 


Initial 
Velocity 
(Vb m/s) 


Initial 

Temperature 
(To K) 


Source 




St-A 


R c = 40 


30 


6.00 


158 


295 


Gust [43] 




St-B 


R c = 40 


30 


6.00 


232 


295 


Gust EH 




St-C 


i? c = 40 


30 


6.00 


183 


715 


Gust EH 




St-D 


i? c = 40 


30 


6.00 


312 


725 


GustE3l 




St-E 


i? c = 40 


30 


6.00 


136 


1285 


Gust (431 




St-F 


R c = 40 


30 


6.00 


160 


1285 


Gust ED 




St-G 


R c = 30 


25.4 


7.62 


208 


298 


Johnson and Cook 


(221 


St-H 


R c = 30 


12.7 


7.62 


282 


298 


Johnson and Cook 


ea 


St-I 


R c = 30 


8.1 


7.62 


343 


298 


Johnson and Cook 


[22] 


St-J 












Addessio et al. (44] 





4.3.1 Room temperature impact: steel 



Figure 11 shows the simulated profile of case St-G without friction. The Johnson-Cook model performs quite well 
in predicting the deformed profile of the specimen. An almost identical profile is obtained if we incorporate friction 
at the impact surface. Similar results are obtained for the other room temperature specimens. 



4 
3.5 

3 
2.5 

2 
1.5 



LX 



Figure 11: Comparison of experimental and computed shape of 4340 steel cylinder (St-G) without friction. The 
axes are in cm. 



4.3.2 High temperature impact: steel 

For high temperature impacts tests, the effect of friction is more obvious in that there is a curling of the edges. 
Figures 12 a),(b),(c),(d) show the simulated profiles of cases St-D and St-F with friction. Specimen D is at a 
lower temperature than specimen F but the impact velocity of the form is almost double that of the latter. The 
Johnson-Cook model predicts the final length of the St-D accurately but underestimates the final length of St-G. 
This indicates that the high temperature behavior of the model is not quite correct even though the rate dependence 
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is captured well. On the other hand, the Steinberg-Guinan model fails miserably at predicting the high velocity 
response but does well for the low velocity/high temperature response. 



Figures 12 e), (f), (g), and (h) show Taylor impact simulations for cases St-D and St-F with particle erosion. 



No significant difference can be seen in the computed profiles when we compare these to the plots in Figure [12} 
except for the SCG model for the St-D sample. Erosion and fracture of the mushroom end does not appear to have 
a first-order effect on the final length of the impact specimen. 
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(c) JC (St-F). 
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(d) SCG (St-F). 
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Figure 12: Comparison of experimental and computed shapes of 4340 steel cylinder with friction. St-D is at 725 K 
and 312 m/s. St-F is at 1285 K and 160 m/s. The axes are in cm. 



5 CONCLUSION 

Lower temperature simulations lead to predicted profiles that are close to those observed in experiment. This is 
true over a range of impact velocities. However, high temperature impacts do not fare so well. 

For the copper specimens, the Johnson-Cook and Steinberg-Guinan models perform better than the MTS model 
both at low and high temperatures. In future work we show that this is partially due to the stress integration 
algorithm used in these calculations. Also, FEM simulations consistently overestimate the final length of specimens 
at high temperatures. 

For the aluminum specimens, the three models, Johnson-Cook, MTS and Steinberg-Guinan, predict accurate 
final lengths and mushroom diameters at room temperature. However, at high temperatures all three models deviate 
from experiment, especially when the strain rate is increased. This indicates a coupling of strain rate and tempera- 
ture that is either not captured by these models or requires further calibration. Mesh dependence due to softening 
appears to be an issue in the MPM simulations. FEM simulations with LS-DYNA predict profiles (at high temper- 
atures) that are less deformed than those predicted by MPM suggesting that the constitutive model evaluation is not 
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as accurate in FEM calculations. 

For the steel specimens, both Johnson-Cook and Steinberg-Guinan perform well at room temperature. How- 
ever, the Steinberg-Guinan model fails under a combination of high impact velocities and high temperatures. The 
Johnson-Cook model is not very accurate at high temperatures but captures rate-dependent effects quite well. Fail- 
ure of the material at the mushroom end does not appear to affect the final length of the specimen significantly. 
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